R Notebook
Quick check of the first no-genetic data demultiplex sequencing from the pilot. NPCs
library(Seurat)The legacy packages maptools, rgdal, and rgeos, underpinning this package
will retire shortly. Please refer to R-spatial evolution reports on
https://r-spatial.org/r/2023/05/15/evolution4.html for details.
This package is now running under evolution status 0
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attaching SeuratObject
Warning message:
R graphics engine version 15 is not supported by this version of RStudio. The Plots tab will be disabled until a newer version of RStudio is installed. library(tidyverse)Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql ── Attaching packages ──────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.4.2 ✔ purrr 1.0.1
✔ tibble 3.2.1 ✔ dplyr 1.1.2
✔ tidyr 1.3.0 ✔ stringr 1.5.0
✔ readr 2.1.3 ✔ forcats 0.5.2
── Conflicts ─────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()#library(DoubletFinder)
library(enrichR)Welcome to enrichR
Checking connection ...
Enrichr ... Connection is Live!
FlyEnrichr ... Connection is available!
WormEnrichr ... Connection is available!
YeastEnrichr ... Connection is available!
FishEnrichr ... Connection is available!
OxEnrichr ... Connection is available!library(clustree)Loading required package: ggraph#library("scClassify")
#library(SingleCellExperiment)
#library("Matrix")
Read in the data
seu <- readRDS("/Users/rhalenathomas/Documents/Data/scRNAseq/ADHD_ZNZ_Mcgill/seu_souporcell_ADHD.rds")
colnames(seu@meta.data) # shows you the meta data
unique(seu$orig.ident)
unique(seu$souporcell_assignment)
unique(seu$souporcell_status)
dim(seu)
table(seu$souporcell_assignment)
Look at the sequencing
VlnPlot(seu, features = c("nCount_RNA","nFeature_RNA"), pt.size = 0.001)
VlnPlot(seu, features = c("nFeature_RNA"), pt.size = 0.001, y.max = 500)
VlnPlot(seu, features = c("nCount_RNA"), pt.size = 0.001, y.max = 50000)
VlnPlot(seu, features = c("nCount_RNA"), pt.size = 0.001, y.max = 5000)
seu$nFeature_RNA %>% summary
seu$nCount_RNA %>% summary
table(seu$souporcell_status)
dim(seu)
Calculate the percent mitochondrial genes
seu <- PercentageFeatureSet(seu, pattern = "^MT-", col.name = "percent.MT")
seu$percent.MT %>% summary
VlnPlot(seu, features = "percent.MT", pt.size = 0.001)
VlnPlot(seu, features = "percent.MT", pt.size = 0.001, y.max = 20)
I’ll use 10% as a cutoff and check if this removes more of the low read cells. I believe this object is already the CellRanger Filtered object.
#Remove any cells with more than 10% mitochondrial counts
#seu.ft <- subset(seu, percent.MT < 10)
#dim(seu)
#dim(seu.ft)
# removed 656 cells
seu.ft <- subset(seu, percent.MT < 20)
dim(seu)
dim(seu.ft)
# removed 339 cells
VlnPlot(seu.ft, features = c("nCount_RNA","nFeature_RNA","percent.MT"), pt.size = 0.001)
VlnPlot(seu, features = c("nCount_RNA"), group.by = "souporcell_status", y.max = 6000)
Warning: Removed 14234 rows containing non-finite values (`stat_ydensity()`).
Warning: Removed 14234 rows containing missing values (`geom_point()`).
VlnPlot(seu, features = c("nFeature_RNA"), group.by = "souporcell_status", y.max = 3000)Warning: Removed 9246 rows containing non-finite values (`stat_ydensity()`).
Warning: Removed 9246 rows containing missing values (`geom_point()`).
I will filter for higher counts and lower counts
# we aren't adding any more filters here
clean up
seu <- seu.ft
rm(seu.ft)
Remove the doublets and the unassigned cells
dim(seu)
[1] 33538 14643Now the workflow to make clusters
seu <- NormalizeData(seu, normalization.method = "LogNormalize", scale.factor = 10000)Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = 2500)Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|var <- VariableFeatures(seu)
top10 <- var[1:10]
p1 <- VariableFeaturePlot(seu)
p2 <- LabelPoints(plot = p1, points = top10, repel = TRUE)When using repel, set xnudge and ynudge to 0 for optimal resultsp2Warning: Transformation introduced infinite values in continuous x-axis
Scale and get PCA
#Linear dimensionality reduction
#Choosing the number of PCs can depend on how many cells you have
seu <- ScaleData(seu)Centering and scaling data matrix
|
| | 0%
|
|=========================== | 33%
|
|===================================================== | 67%
|
|================================================================================| 100%seu <- RunPCA(seu, assay = "RNA", npcs = 30)PC_ 1
Positive: CKB, TUBB2B, CRABP1, TUBA1A, MAP1B, TTYH1, SOX2, POU3F2, BEX1, CRMP1
HMGA1, KIF5C, MAP2, RFX4, HES6, SOX11, DRAXIN, FGFBP3, PRTG, SRGAP3
HSP90AA1, TAGLN3, MIAT, LIN28A, GPM6B, SYT11, MAP6, ELAVL2, DCX, ELAVL3
Negative: COL3A1, S100A11, LGALS1, LUM, DCN, SPARC, COL1A1, PCOLCE, B2M, CDH11
SERPINH1, APOE, COL5A2, FN1, COL5A1, COL6A2, IFITM3, RRBP1, TIMP1, SELENOM
COLEC12, BGN, COL6A3, LAPTM4A, CALD1, COL1A2, EMP3, IFI16, TPM1, CCDC80
PC_ 2
Positive: HMGB2, CENPF, TOP2A, CKS1B, MKI67, NUSAP1, BIRC5, TPX2, HMGB1, CCNB1
H2AFZ, ASPM, TUBA1B, H2AFX, UBE2C, CKS2, SMC4, MAD2L1, DLGAP5, GTSE1
CENPE, PTTG1, TYMS, TUBB4B, PCLAF, CDK1, CCNA2, CDC20, KIF20B, KIF11
Negative: STMN2, DCX, TAGLN3, ELAVL3, INA, STMN4, NEFM, KLHL35, ELAVL4, NHLH1
ONECUT2, SCG3, DCC, SRRM4, SYT1, INSM1, DLL3, CNTN2, MLLT11, CDKN1C
NCAM1, TFAP2B, CLDN5, ATCAY, SSTR2, SOX4, GDAP1L1, GNG3, ELAVL2, PCBP4
PC_ 3
Positive: IFI44L, PRRX1, SMOC2, PTH1R, COLEC12, SFRP1, LRRC17, CDC42EP5, DAB2, TWIST1
SLC1A3, TPM1, ZIC1, CPE, NAV3, FOXC1, IFI44, FRZB, COL26A1, ID1
DACT2, PRRX2, LAMA4, FN1, TCIM, PMP22, EDNRA, ISG15, AKAP12, NFIA
Negative: KRT19, KRT8, THY1, TM4SF1, NPNT, S100A10, KRT18, HAND2, MYRF, NNMT
DSG2, MATN2, GATA6, FXYD5, RPP25, ASPHD1, PLP2, SLC9A3R1, KCNK6, LRRN4
PITX1, MAN1A1, CYBA, LMCD1, KDR, EZR, TAGLN2, COL1A1, BNC1, ANXA1
PC_ 4
Positive: NEB, DES, KLHL41, TNNI1, TNNT1, TTN, CHRNA1, MYL1, MYL4, MYOG
TNNC2, ENO3, CDH15, MYOD1, ACTC1, TNNT2, MYLPF, SMPX, SGCA, AC020909.2
SYNPO2L, IL17B, MRLN, FITM1, MYH3, ATP2A1, SMYD1, TRIM55, TNNC1, RAPSN
Negative: PTN, IGFBP2, CRABP1, IFI16, GJA1, TTYH1, PDLIM1, ZIC2, COL11A1, CLU
SOX2, FABP5, IFIT3, MGP, PARP14, SFRP2, IFITM1, SERPING1, VCAN, SP100
HLA-C, SFRP1, C7, B2M, DAB2, ARHGAP29, IFI6, HLA-E, GPM6B, LIFR
PC_ 5
Positive: RPL12, VIM, KRT19, MYLPF, IGFBP5, MYL1, LMCD1, TNNC2, AL359091.1, MYL4
HAND2, HES5, TNNC1, NPNT, NPPC, FXYD5, PITX1, NUPR1, KRT8, ACTC1
TNNI1, TNNT2, MYRF, SPRR2F, AL139246.5, IL17B, TTYH1, HSPB7, REC8, RPP25
Negative: ASPM, ARL6IP1, CENPE, NEFM, TPX2, UBE2S, GAP43, STMN2, TOP2A, MKI67
KPNA2, NCAM1, INA, GTSE1, SGO2, UBE2C, STMN4, ELAVL3, MAP1B, DLGAP5
NEFL, NUSAP1, PLK1, CENPF, ONECUT2, CCNB1, DCX, ACTB, KIF20B, HMMR PCAPlot(seu)
#Assess how many PCs capture most of the information in the data
ElbowPlot(seu, ndims = 30)
20 PCs are good now make the UMAP
seu <- RunUMAP(seu, dims = 1:20, n.neighbors = 121)Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
15:34:51 UMAP embedding parameters a = 0.9922 b = 1.112
15:34:51 Read 14643 rows and found 20 numeric columns
15:34:51 Using Annoy for neighbor search, n_neighbors = 121
15:34:51 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:34:53 Writing NN index file to temp file /var/folders/k4/khtkczkd5tn732ftjpwgtr240000gn/T//Rtmpo5qc3e/file10214368d353
15:34:53 Searching Annoy index using 1 thread, search_k = 12100
15:35:09 Annoy recall = 100%
15:35:09 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 121
15:35:13 Initializing from normalized Laplacian + noise (using irlba)
15:35:14 Commencing optimization for 200 epochs, with 1862028 positive edges
Using method 'umap'
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:35:36 Optimization finishedDimPlot(seu, group.by = "orig.ident")
Find clusters
seu <- FindNeighbors(seu, dims = 1:20, k.param = 121)Computing nearest neighbor graph
Computing SNN# the number of clusters is dependent on the resolution a number from 0-2.
# Higher values make more clusters
# we include
seu <- FindClusters(seu, resolution = c(0,0.05,0.25,0.6,1,1.5) )Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 14643
Number of edges: 2834426
Running Louvain algorithm...0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|Maximum modularity in 10 random starts: 1.0000
Number of communities: 1
Elapsed time: 8 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 14643
Number of edges: 2834426
Running Louvain algorithm...0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|Maximum modularity in 10 random starts: 0.9746
Number of communities: 5
Elapsed time: 12 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 14643
Number of edges: 2834426
Running Louvain algorithm...0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|Maximum modularity in 10 random starts: 0.9249
Number of communities: 7
Elapsed time: 9 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 14643
Number of edges: 2834426
Running Louvain algorithm...0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|Maximum modularity in 10 random starts: 0.8785
Number of communities: 10
Elapsed time: 13 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 14643
Number of edges: 2834426
Running Louvain algorithm...0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|Maximum modularity in 10 random starts: 0.8384
Number of communities: 13
Elapsed time: 11 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 14643
Number of edges: 2834426
Running Louvain algorithm...0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|Maximum modularity in 10 random starts: 0.8012
Number of communities: 20
Elapsed time: 9 seconds# we can visualize which cells are grouped together at different resolutions using clustree
clustree(seu, prefix = "RNA_snn_res.")
# 0.6 looks good to annotate. Each cluster is splitting apart up to this point and then the cells start merging and changing clusters.
Visualize
DimPlot(seu, group.by = "RNA_snn_res.0.6", label = TRUE)
DimPlot(seu, group.by = "RNA_snn_res.1", label = TRUE)
# clustering could be improved
Look at expression profiles of known cell type markers
da_neurons <- c("TH","SLC6A3","SLC18A2","SOX6","NDNF","SNCG","ALDH1A1","CALB1","TACR2","SLC17A6","SLC32A1","OTX2","GRP","LPL","CCK","VIP")
NPC_orStemLike <- c("DCX","NEUROD1","TBR1","PCNA","MKI67","SOX2","NES","PAX6","MASH1")
mature_neurons = c("RBFOX3","SYP","DLG45","VAMP1","VAMP2","TUBB3","SYT1","BSN","HOMER1","SLC17A6")
excitatory_neurons = c("GRIA2","GRIA1","GRIA4","GRIN1","GRIN2B","GRIN2A","GRIN3A","GRIN3","GRIP1","CAMK2A")
inhbitory_neurons = inh = c("GAD1","GAD2", "GAT1","PVALB","GABR2","GABR1","GBRR1","GABRB2","GABRB1","GABRB3","GABRA6","GABRA1","GABRA4","TRAK2")
astrocytes <- c("GFAP","S100B","AQP4","APOE", "SOX9","SLC1A3")
oligodendrocytes <- c("MBP","MOG","OLIG1","OLIG2","SOX10")
opc <-
radial_glia <- c("PTPRC","AIF1","ADGRE1", "VIM", "TNC","PTPRZ1","FAM107A","HOPX","LIFR",
"ITGB5","IL6ST","SLC1A3")
epithelial <- c("HES1","HES5","SOX2","SOX10","NES","CDH1","NOTCH1")
microglia <- c("IBA1","P2RY12","P2RY13","TREM119", "GPR34","SIGLECH","TREM2",
"CX3CR1","FCRLS","OLFML3","HEXB","TGFBR1", "SALL1","MERTK",
"PROS1")
features_list <- c("MKI67","SOX2","POU5F1","DLX2","PAX6","SOX9","HES1","NES","RBFOX3","MAP2","NCAM1","CD24","GRIA2","GRIN2B","GABBR1","GAD1","GAD2","GABRA1","GABRB2","TH","ALDH1A1","LMX1B","NR4A2","CORIN","CALB1","KCNJ6","CXCR4","ITGA6","SLC1A3","CD44","AQP4","S100B", "PDGFRA","OLIG2","MBP","CLDN11","VIM","VCAM1")
short_list <- c("MKI67","SOX9","HES1","NES","DLX2","RBFOX3","MAP2","TH","CALB1","KCNJ6","SLC1A3","CD44","AQP4","S100B","OLIG2","MBP","VIM")
View on UMAP
Idents(seu) <- "RNA_snn_res.0.6"
for (i in NPC_c) {
print(FeaturePlot(seu, features = i, min.cutoff = 'q1', max.cutoff = 'q97', label = TRUE))
}Warning in FeaturePlot(seu, features = i, min.cutoff = "q1", max.cutoff = "q97", :
All cells have the same value (0) of FOXG1.
Idents(seu) <- "RNA_snn_res.0.6"
for (i in astrocytes) {
print(FeaturePlot(seu, features = i, min.cutoff = 'q1', max.cutoff = 'q97', label = TRUE))
}
